2002 - 01-2333 


Towards a 3D Space Radiation Transport Code 

J.W. Wilson, R.K.Tripathi 

NASA Langley Research Center, Hampton, VA 

F.A. Cucinotta 

NASA Johnson Space Center, Houston, TX 

J.H. Heinbockel, J. Tweed 

Old Dominion University, Norfolk, VA. 


Copyright© 2002 Society of Automotive Engineers, Inc. 


ABSTRACT 

High-speed computational procedures for space 
radiation shielding have relied on asymptotic expansions 
in terms of the off-axis scatter and replacement of the 
general geometry problem by a collection of flat plates. 
This type of solution was derived for application to 
human rated systems in which the radius of the shielded 
volume is large compared to the off-axis diffusion limiting 
leakage at lateral boundaries. Over the decades these 
computational codes are relatively complete and lateral 
diffusion effects are now being added. The analysis for 
developing a practical full 3D space shielding code is 
presented. 

INTRODUCTION 

Early methods of space radiation shield evaluation relied 
largely on Monte Carlo codes 11,2] and made important 
contributions to NASA engineering programs. Yet, slow 
computational procedures did not allow progress in 
coupling to spacecraft geometry for over 30 years 
(compare Alsmiller et ai. [3], Armstrong et al. [4]) and did 
not allow early entry of radiation constraints into the 
design process and off-optimum Monte Carlo solutions 
to shielding problems continue to plague final designs 
[5-7]. Simulations with full 3D Monte Carlo codes often 
use questionably simplified ID shielding geometry 
models to increase computational speed in which 
shielding for an ISS module is approximated as a very 
long 20.7 g/cm 2 thick aluminum cylindrical shell [4] 
leading to an overestimate of the neutron flux within ISS 
since neutron leakage is underestimated in these 
homogenized configurations. The "3D calculations” 
used in current ISS Monte Carlo studies are routinely 
exceeded in their correctness in coupling to spacecraft 
geometry by HZETRN [6-8]. In addition to spacecraft 
geometry simplifications, the astronauts in Monte Carlo 
calculations are approximated by a 30-cm sphere of 
tissue [3,4]. Conversely, the use of a 3D Monte Carlo 


code within such a simplified ID geometry for ISS will 
lead to erroneous underestimates of the trapped 
charged particle exposure since anisotropic shield 
distributions within ISS have long been known to be a 
major factor (2 to 5) in astronaut organ exposure [9], 

The reason for such simplifications using Monte Carlo 
methods is easy to understand. To map the interior of an 
ISS module requires at least 10 6 events within each 
volume about a field point. With 10 3 points for a field 
map, this requires 10 9 events to map an interior field. 
Each particle track from boundary to exit will experience 
1 0 3 -1 O'* events requiring a ray trace per event leading to 
10 12 -10 13 ray traces for a field evaluation. In our 16A ISS 
configuration in which only the HAB module is currently 
represented in detail [6], 14 seconds is required to 
evaluate 10 3 rays with a 12 processor script on an Origin 
parallel machine or approximately 14 msec per ray using 
optimized ray tracing methods [10]. The Monte Carlo 
interior field evaluation from ray tracing alone on a 1 2 
processor Origin would require lO’MO 11 seconds. It 
should be clear that Monte Carlo methods would not be 
useful in design optimization. 

The development of high-speed deterministic 
procedures allows early entry of radiation constraints into 
the design process without simplification of the shield 
geometry [1 1] but Monte Carlo methods could still play a 
role in future final design evaluation with full geometry 
provided more efficient ray trace methods can be found, 
massively parallel machines are used, and/or a radically 
different approach to ray tracing can be found. Current 
Monte Carlo methods simply use trivial geometry, which 
is computationally efficient [4], but provide generally 
incorrect results and do not meet the requirement for 
efficient linkage to spacecraft geometry or biological 
models. 

In the present report, we will give an overview of the 
salient features allowing the development of high- 


performance computational procedures for evaluation of 
the 3D geometry of spacecraft with astronauts. All of the 
steps in the process have been demonstrated at some 
minimal conceptual level and much of the software is 
available in terms of a computationally efficient first order 
approximation. Only the final assembly of the software 
package with the spacecraft/astronaut geometry and 
further evaluation of the cross section database is 
required. This last step will require a well-defined 
laboratory and flight validation of the program. 
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Fig.1. Schematic of a point source of monoenergetic 
partilcles 
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Fig. 3. Profile of 400 MeV neutrons. 


PRELIMINARY CONSIDERATIONS 

Early Monte Carlo studies revealed solutions to nucleon 
transport that showed dose to be confined near the 
particle path as shown in fig. 1. The placement of near 
boundaries leads to leakage but is ineffective except 
near the central region as for boundary B. Further 
studies using an early version of HETC [12] reveals 
similar results for high-energy proton and neutron beams 
as seen in figs. 2 and 3. That even high-energy neutron 
beams are confined to a few centimeters of the beam 
axis results from the dominance of low-momentum 
transfer charge-exchange processes. Outside the few 
centimeters, the dose drops by an order of magnitude 
resulting from the diffuse low-energy neutron 
component. 

Other testing with HETC for broad beams of normal and 
isotropic incidence protons showed only few percent 
differences between a straightahead approximation and 
full 3D calculations [11,13] and provided a starting point 
for development of high-speed computational methods 
provided in part by the rapid convergence of the 
Neuman series for proton transport that tends to be 
dominated by atomic interactions at medium energies 
[14]. The neutron/Neuman series converges more 
slowly due to the neutron charge neutrality [14], One 
could show that ID errors committed in 3D-space 
applications was on the order of the second power of the 
ratio of beam divergence to the radius of curvature of 
space shields [11] which produces little error for human 
rated systems. 

Adding 3D multiple-elastic scattering to the Neuman 
series allowed good comparisons with time-of-flight 
experiments with 670 A MeV Ne beams where lateral 
diffusion dominated the measurements (few percent 
detector acceptance [15,16]. Further laboratory 
validation with broad C, N, O, and Fe ion beams has been 
completed giving confidence to our understanding of 
ion interactions in bulk materials [7,16-23]. The first-order 
Neuman series provided the basis for a high-speed 
marching procedure for coupled proton/neutron/HZE 
transport (HZETRN) in common usage in space radiation 
protection [24], Target fragment contributions are 
evaluated analytically within the computational procedure 
[24-29], 

Only the lateral scattering at the fragmentation event and 
the lateral diffusion of light ions and especially neutrons 
remains to be evaluated. Several approaches have been 
followed for a full 3D treatment of coupled neutron 
transport, a standard Sn code with a multigroup cross 
section set from reactor applications using first order 
finite difference methods was employed by coupling to 
HZETRN [30] and a separate line of development first 
examining mean-value multigroup methods without finite 
differencing and methods of collocation [31]. The 
HZETRN/neutron multigroup code has had reasonable 
success in shuttle flight validation [8]. A parallel version 



of HZETRN which operating on a 10 processor Alpha 
chip cluster is able to solve the Boltzmann equation for 
the HAB module crew quarter augmentation in few 1 0’s 
of minutes which is sufficient for implementation of 
optimization procedures now being developed [6]. 

At this juncture, the deterministic evaluation of multiple- 
elastic scattering of heavy ions has been adequately 
resolved with detailed laboratory testing [15,16,32]. The 
perplexity of straggling has been resolved [33]. The 3D 
geometry in charged particle transport can be accurately 
approximated by added asymptotic terms [1 1 ,34]. The 
lowest order asymptotic term have been solved as a 
marching procedure [24], a perturbation expansion [18], 
and a nonperturbative expansion [35], with extensive 
testing in laboratory and flight experiments. Detailed 
investigations of the diffusive neutron component at 
various levels of approximation from 3D Sn methods to 
collocation methods including PN methods have shown 
promise in approaching an efficient full 3D code but 
needs further development. It has been demonstrated 
that these methods can be used efficiently in the full 3D 
complexity of spacecraft design and have been flight- 
tested. Added treatment of the charged particle 
diffusive components, which will converge more rapidly 
than the neutral neutrons can be adequately treated 
using perturbation theory. A meson database needs 
completion [36] with addition of a complete E&M 
cascade code [37]. Further considerations of fast 
computational procedures and implementation issues on 
high-performance computers for ultimate insertion into 
shield optimization and reliability design methods needs 
to be made. 

3D CODE DEVELOPMENT 

The relevant transport equations are the linear 
Boltzmann equations for the flux density 0j(x, £2, E) for 
particle type j as 

&v<t> j(x,o,E)- 2fofr(Q,ar,E,E') Mx,ar,E) dE' 

- oj(E) *j(x,fl,E) (1) 

where oj(E) and ojk(G,i2’,E,E') are the shield media 
macroscopic cross sections. The ojk^C.E.E') 
represent all those processes by which type k particles 
moving in direction £2' with energy E’ produce a type j 
particle in direction £2 with energy E (including decay 
processes). Note that there may be several reactions 
that produce a particular product, and the appropriate 
cross sections for equation (1) are the inclusive ones. 
Exclusive processes are functions of the particle fields 
and may be included once the particle fields are known. 
The total cross section cj (E) with the medium for each 
particle type is 


where the first term refers to collision with atomic 
electrons, the second term is for elastic nuclear 
scattering, and the third term describes nuclear reactions 
where we have ignored the minor nuclear excitation 
processes. The corresponding differential cross section 
is similarly ordered. Many atomic collisions (~ 10®) occur 
in a centimeter of ordinary matter, whereas - 10 3 nuclear 
coulomb elastic collisions occur per centimeter, while 
nuclear reactions are separated by a fraction to many 
centimeters depending on energy and particle type. 
Solutions use the atomic collisions as first perturbation 
with special methods used for neutrons for which atomic 
cross-sections are taken as zero. 

The cross sections appearing in the Boltzmann equation 
are the inclusive ones so that the time-independent 
fields contain no spatial (or time) correlations and are 
evaluated once the fields are known [24]. Such 
correlations are important to the biological injury of living 
tissues [38] and transients in devices [27]. For example, 
the correlated release of electrons or target fragments 
about the ion's path determines the biological response 
to these ions [26,28,38] and is the driver of shield 
material attenuation properties [39]. The solution of 
equation (1) involves hundreds of multi-dimensional 
integro-differential equations that are coupled together 
by thousands of cross terms and must be solved self- 
consistently subject to boundary conditions. 

The first order physical perturbation to the right side of 
equation (1) is the atomic/molecular cross sections as 
noted in equation (2) for which those terms in (1) are 
expanded about the energy moments resulting in 
range/energy relations including relativistic polarization 
effects [40] and straggling parameters [21,24,33]. 

The distribution of the electrons about the ion path is 
critical to evaluation of biological injury [38], is critical to 
the evaluation of shield attenuation properties [39], and 
fundamental to dosimetric evaluation of astronaut 
exposure risks [41]. Such effects are likewise governed 
by equation (1). The next physical perturbation term is 
the coulomb scattering by the atomic nucleus and is 
represented by Rutherford scattering modified by 
screening of the nuclear charge by the orbital electrons 
using the Thomas-Fermi distribution for the atomic 
orbitals. The total nuclear coulomb cross section found 
by integrating over the scattering directions is related to 
the radiation length [32] and has been well validated in 
laboratory tests [15,16,32]. The electron production in 
ion collisions and the electromagnetic cascade following 
the neutral meson decay are described, in part, 
elsewhere [36,37]. 

The nuclear reactive cross sections can be written in the 
following form 


(E) = °j,at (E) + cj.el (E) + c ), r (E) 


( 2 ) 


a^ r (£2,£2 ,E,E ) - c^k,» 0 (E,E ) + Op. far (£2,£2 ,E,E ) 


( 3 ) 


where the first term is isotropic and associated with lower 
energy particles produced including target fragments 
and the second term is highly peaked in the forward 
direction and is associated mainly with direct quasi-elastic 
events, charge exchange, and projectile fragments. 
Surprisingly even nucleon-induced reactions follow this 
simple form and the isotropic term extends to relatively 
high energies (see Fig. 4). For nucleon induced 
reactions, the following form has been used in versions 
of FLUKA as follows 

<kAQ#JE,E) = v jK (E’) <v(E’)f jk (E.E’) g R (Q*a,E,A r ) (4) 

where the Ranft factor used in FLUKA is 

q h (£ 2*G’,E,A t ) - N r expI-^/AJ nfodzO (5) 

and constant for larger values of production angle 0and 
given by Ranft as 

X R = (0.12 + 0.00036A t /E) ( 6) 

although new generalized fits are being derived. This 
separation in phase space is exploited for efficient 
computational procedures. The heavy-ion projectile 
fragment cross-sections are further represented by 

o^Q,a, E,E') = oj kr (E’) N, exp[-2mvtEE ’)(l-Q'0)le IJk ] 
xexp[ - (E + \-E') 2 /2 tf/S(2xEj) (7) 

where A, k is related to the momentum downshift, % is 
related to the longitudinal momentum width, e IJk is 
related to the transverse momentum width, and N, is the 
transverse normalizing factor. Since the transverse width 
is small compared to the projectile and fragment energy 
the transverse function is highly peaked about the 
forward direction (ie, - 1). 



Atomic interactions limit the contributions of charged 
particles in the transport process. For example, the 
protons and alpha particles produced in aluminum below 
100 A MeV contribute to the fluence only within a few 
centimeters of their collision source and the heavier ions 
are even more restricted (see Fig. 5). This is an important 
factor in that the transported secondary charged particle 
flux tends to be small at low energies and the role of 
additional nuclear reactions are likewise limited. 



vector array field function as 

4>=fo(x.aE)] (8) 

the drift operator 

D = [Q*V\ (9) 


Fig. 4. Isotropic and forward neutron spectra 
produced bv 500 MeV proton in aluminum. 


and the interaction operator 


I = \Ifa^{Q,Sf,E,E) di? dE' - 05(E)] 


( 10 ) 




with the understanding that I has three parts associated 
with atomic, elastic, and reactive processes as given in 
equation (2). Equation (1) is then rewritten as 

= ■,•* (11) 

where the first two physical perturbation terms are shown 
on the left-hand side and have been adequately 
resolved in past research. The reaction cross section is 
separated by equation (3) into isotropic and forward 
component for which equation (11) may be written as 
coupled equations 

[D-I. t — l.i + <*>«« = 

E,E') d£2* dE’}* 0 for - E rM <P tor (12) 

and 

[D - I* -l,i+ aj* *,..={/ o r (£2,£2*, E,E') d Q' dE’}* *,. e 
+ {/V J£2,£2\ E,E*) d£2* dE'}» <P tor 

* “r * ^t .0 + “r,i«o * ^tor (13) 

Equation (12) can be solved as a Neuman series [24,34] 
and written as 

^for = [G + G , -r,for*G + G *“r,»or *G *“r,for *G + (14) 

for which the series can be either evaluated directly or 
proscribed as a marching procedure in either a 
perturbative sense as the current form of HZETRN or 
nonperturbative sense (future version of HZETRN) as 
described elsewhere [42], 

The cross term in equation (13) gives rise to an isotropic 
source of light ions and neutrons of only modest 
energies for which Fig. 4 is typical. Spectral contributions 
to the Neuman series depend on the particle range and 
probability of surviving nuclear reactions that establish 
the functional form of the G matrix. The second term of 
the Neuman series is proportional to the probability of 
nuclear reaction that is limited by the particle range as 
discussed above and shown in Fig. 6. It is clear from Fig. 
6 that nuclear reactions for the charged particles below a 
few hundred A MeV are infrequent for which fast 
convergence has been demonstrated. For the moment 
we will neglect the straggling and multiple-elastic 
processes to simplify the present explanation (provide 
only minor corrections to space radiation exposures) and 
examine the remaining reactive terms of equation (12). 
The corresponding Volterra equation is given by 

<P\M E) = (Sj(E r )Pj(E r ) <Pi(r(Q,x)M E r ) 

+ ^dE’A^EUXdEW o jkfor (£2,& ,E’,E”) 
x ^(x+lR^E) - Rj(E')]O t i7 ,E”)}/ Sj(E)Pj(E) (15) 

where r is the point on the boundary connected to x 
along -£2, E, = Rj -1 [p - d + RJ, p is the projection of x onto 
£2, and d is the projection of r onto £2. Equation (14) 


results from the Neuman series solution to equation (15). 
In the past we have expanded the angular integral Q' 
asymptotically and implemented as a marching 
procedure (HZETRN, [43]), as a perturbation expansion 
[18], and by non-perturbative approximation [44] 
resulting in three distinct methods to evaluate the first 
order asymptotic terms, all of which have had extensive 
experimental validation. Independent of the method 
used to evaluate the lowest asymptotic order term, the 
first correction term is found by replacing the fluence in 
the integrand of equation (15) by the lowest order 
asymptotic solution as 

*j(x,fl,E) = {Sj(E..)P,(E.) <pi(r(Q,x),Q,E r ) 

+ ^dE'A^EyX dEW o jk JG,a,E’,E") 
x 0 ko (x+[Rj(E) - R^EW.^E”)}/ Sj(E)Pj(E) (16) 

where 0j(x,£2,E) is found as an integral over the 
neighborhood of rays centered on £2 using the lowest 
order asymptotic solution 0 ko (x,£2’,E”) along an adjacent 
ray directed along £2’. Note that the boundary condition 
reached along -£2’ enters through the lowest order 
asymptotic approximation and the angular integral 
correction in equation (16) is determined by the 
homogeneity and angular dependence of the space 
radiation and radius of curvature of the bounding material 
which we have shown long ago are the determinant 
factors of the magnitude of the first order asymptotic 
correction which is anticipated to be very small for human 
rated systems (large radius of curvature) in space 
radiation which is homogeneous and isotropic in most 
applications [11,24,34]. 

In a region of small radius of curvature the specific flux 
components near the site of evaluation will be missing 
contributions along adjacent rays which do not 
compensate losses along the ray on which the solution is 
evaluated representing the losses due to leakage. This 
computational procedure is only a small addition to prior 
code development and will have little impact on 
computational efficiency. The angular dependence of 
the integral kernel of equation (16) is controlled by the 
forward reactive cross section o^^C.E’.E’’) with 
its highly peaked structure given by equations (4) or (7) 
depending on particle type. Even the most dispersed 
(nucleonic) components remain near the forward 
direction as shown in Figs. 2 and 3. The angular 
dependence of the forward peak of fragmenting Ca ions 
at 100 and 1,000 A MeV is shown in Fig. 7. The low- 
energy ions with limited range have transverse 
components on the order of 10 degrees reducing to a 
few degrees at high energies. Note that the low energy 
ions have limited range and will contribute little to the 
transported flux (see Fig. 5) or nuclear reactions (see Fig. 
6). The higher energy ions with their much longer path- 
lengths giving more important contributions are related 
to only a very small angle of acceptance (few degrees) at 
the boundary. The form of the kernel leads directly to a 



* Gauss-Hermite expansion and evaluation over the angle 
of production. Although the neutron Neuman series for 
the forward components converge more slowly since 
their contribution to the neutron flux is not limited by 
atomic interactions these higher energy neutrons will be 
adequately evaluated by the same procedures (for 
example, see Fig. 3). Higher order asymptotic terms can 
be evaluated with similar iteration of equation (15) if 
required but all indications are that the first such 
correction will be small [1 1 ,34], 
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Fig. 7. Normalized transverse components for Ca 
fragmentation. 

The remaining problem is solution for the transport of the 
low-energy neutron and light ion isotropic sources in 
equation (1 3) that dominate the solution below about 70 
A MeV (see Fig. 4). In this region light ion transport is 
completely dominated by the atomic interaction terms 
and only a very small fraction have nuclear reactions 
making only minor contributions to the particle fields (Fig. 
6). This is especially true for the target fragments that 
can be solved in closed form [25,34], The 3D solution 
for the low-energy neutrons and light ions is approached 
through the integral of equation (16) with the inbound 
flux equal to zero as 

f|(x,fl, E) = 

OT'dE'AfitEU"/* dEW a Jk jai?,E’,E”) 

* k (x+[Rj(E) - Rj(E’)]i2,C ,E”) +V E Er dE’A,P j (E*) 
S k ( x+ [Rj(E) - Rj(E’)]£?,E')}/ Sj(E)Pj(E) (17) 

where ^(x,E’) is the isotropic source from coupling to 
equation (13). We solve equation (17) by perturbation 
where the lowest order term is only the source term. We 
evaluate the first correction by replacing 0 k (x+[Rj(E) - 
under the integral by the lowest order 
solution. The correction is proportional to the nuclear 


reaction probability as given in Fig. 6 and is quite small 
except for neutrons. 

The neutrons have no charge and are 
dominated at low energies by elastic and reactive nuclear 
processes. We further expand equation (13) for the 
single neutron component as 

[Q'V + oj 0n(x,aE) -/CTnn(£?,tf,E,E') 0n(x,i7,E’) 
d£? dE' + [4*o» <fw+ 4* ^(ujn (18) 

where the last term is from coupling to the solution of 
equation (17) and the spectral properties are 
characterized by 

f,(E) =/f(E,E’) UE’,500M eV) dE’ (19) 

with results shown in Fig. 8. The neutron spectrum is 
greatly degraded in energy on the first collision (Fig. 8) 
and the remaining development of the low-energy 
neutron transport is the last issue to be resolved. This is 
the typical nuclear engineering problem for which a 
multitude of methods have been developed such as the 
Sn, multigroup, and collocation methods already applied 
to HZETRN. It is mainly a question of computational 
efficiency and we continue to investigate this issue. 



Figure. 8. Second collision neutron spectrum fl(E). 


Although, equation (18) exhibits behavior similar to 
thermal diffusion there are strong differences between 
thermal and neutron diffusive processes. Thermal 
diffusion at ordinary temperatures has minor leakage 
through near boundaries since radiative processes are 
proportional to T 4 (in the absence of convection) leaving 
lateral diffusion an important process. In distinction, 
neutron diffusion is dominated by leakage at near 
forward and backward boundaries in penetrating 
spacecraft walls and lateral diffusion play a minor role. 
Generally, low-energy neutron leakage is a dominant 
process within 15-20 g/cm 2 of the bounding surface in 




* most materials. Since human rated systems have 
shielding of large radius of curvature and small thickness 
to radius ratio as determined by living and working space 
requirements, it approximates a connected system of flat 
plates for which leakage at forward and backward 
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Fig. 8. Observed and calculated proton spectra. 
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Fig. 8. Observed and calculated deuteron spectra. 



Fig 10. Observed and calculated neutron spectra. 


boundaries dominates the transport. In this limit the 
transport simplifies to a connected set of ID transport 
problems with leakage at the back and forward 
boundaries [8,31] and shows reasonable success in 
comparison with experimental flight data, see Figs. 8-10 
[45,46]. In the present development we considered a 
convergent series of approximations to gauge accuracy 
of the transport procedures and allow choices of the 
most practical method in a given application. The results 
shown in Figs. 8-10 are the lowest order asymptotic term 
and the good comparison with flight data encourages us 
to believe that the earlier analysis on the rates of 
convergence of the asymptotic series [11,34,35] is 
correct and the next corrections will be small with the 
possible exception of improved low energy neutron 
transport results. We will of course require improved 
neutron measurements as seen from the two analyses of 
the JSC bonner sphere data. 

CONCLUSION 

Considerable research on the development of high- 
performance computational procedures has been 
committed over the last thirty years. The problem of 
developing a full 3D version of these procedures is now 
very clear and the steps required have been fully 
demonstrated. Even a partial solution of the problem has 
met with considerable success in comparison with 
laboratory and flight measurements. A program for 
developing the final software is now under way and 
should be complete over the next several years. 
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